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We sampled 45 Andean mountain tapirs (Tapirus pinchaque) from Colombia and Ecuador and sequenced 
15 mitochondrial genes (two rRNA and 13 protein codifying genes)—making up 13,939 base pairs, approx¬ 
imately 83.1% of the total mitochondrial DNA’s length. The overall sample had low to medium levels 
of nucleotide diversity with diversity slightly higher for the Colombian population. Both populations 
experienced high historical gene flow and our genetic heterogeneity analyses revealed a low genetic 
differentiation between them. Therefore, we did not detect any molecular subspecies, or significantly 
different evolutionary units for T. pinchaque. This species experienced a population expansion in the last 
100,000 years but this expansion was more pronounced in the Ecuadorian population especially in the 
last 10,000 years, whereas the Colombian population underwent a strong bottleneck in the last 5,000 
years. There was no significant spatial trend in genetic structure for the mountain tapir in Colombia 
and Ecuador. Phylogenetic analyses did not detect any important geographic clade within this species. 
Temporal split between T. pinchaque and T. terrestris might have occurred around 7-1.5 million years 
ago (MYA). T. pinchaque and T. terrestris + T. kabomani are two monophyletic clades, suggesting that T. 
kabomani is not a full species. 

© 2015 Deutsche Gesellschaft fur Saugetierkunde. Published by Elsevier GmbH. All rights reserved. 


Introduction 

An important first step in protecting groups of similarly struc¬ 
tured organisms is to place them into discrete taxa. Once a species 
is recognized we monitor its overall fitness and develop conser¬ 
vation plans towards its benefit. Unfortunately, biologists do not 
always agree on the naming of new species in part due to the 
species concept they select (i.e. Phylogenetic Species Concept, PSC, 
vs. Biological Species Concept, BSC) as well as the data (fossils, DNA 
results, morphology, ethology, ecology) they use. One example of 
this non-agreement is in the genus Tapirus. In Latin America three 
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tapir species have been traditionally accepted, the Baird or Central 
American tapir (Tapirus bairdii ), the mountain tapir (T. pinchaque ) 
and the lowland or Brazilian tapir (T. terrestris). There is also the 
Malayan tapir species (T. indicus) in Asia (Brooks et al., 1997). 

T. pinchaque (Roulin 1829) is the smallest of the four living tapir 
species—having an average shoulder height of 0.90 m, body length 
of 1.8 m and weight of 150 kg. This species is found in habitats of 
1,700-4,800 m above sea level (masl) (Downer, 2001; Schauenberg, 
1969) in different types of north Andean forests and paramos (open 
grass lands) in Colombia, Ecuador, and in Northern Peru. It is an 
extremely efficient seed disperser. Downer (2001) documented 
more than 50 Andean plant species that have germinated in tapir 
feces within the Sangay National Park (Ecuador). Mountain tapir 
survival is considered a crucial factor for the conservation of the 
northern Andean wilderness and watershed, an area containing the 
most important water reserves in this and surrounding areas. 


http://dx.doi.Org/10.1016/j.mambio.2015.ll.001 
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This tapir is currently classified as an endangered species in the 
IUCN Red List (1UCN 2011) at the global level and it is also enclosed 
in Appendix I in CITES, being one of the rarest large mammals of the 
world. Moreover, they have been categorized as critically endan¬ 
gered in the Andean area of Colombia, Ecuador and Peru (Lizcano 
et al., 2006; Tirira, 2011 ). 

The origin of tapirs in South America and the role played by T. 
pinchaque in this origin are very controversial topics. Hershkovitz 
(1954) analyzed morphological characters and concluded that the 
ancestor of T. pinchaque was the first tapir form to enter South 
America prior to the formation of the Panamanian Isthmus more 
than 3-4 MYA. Later, the ancestors of T. terrestris, T. bairdii and other 
tapirs from hotter climates migrated into South America. Therefore, 
the ancestors of the three tapir species currently present in South 
America previously appeared in North America at a time when the 
Isthmus of Panama was not yet formed. Three independent lineages 
evolved and, according to Hershkovitz (1966), T. bairdii is the most 
recent. This author distinguished two subspecies of T. pinchaque: T. 
p. pinchaque (Roulin; type locality in Paramo de Sumapaz in Cun- 
dinamarca, Colombia) and T. p. leucogenys (Gray; type locality in 
Paramo de Azuay in the Eastern Cordillera of Ecuador). 

Haffer (1970) introduced a different hypothesis on the origin 
of the tapirs of South America taking into account morphological 
and biogeographical considerations. He considered the ancestor 
of T. pinchaque to be the first form which colonized the North¬ 
ern Andes during its last rising within the Pliocene. During the 
Pleistocene, the ancestor of T. pinchaque gave rise in the eastern 
lowlands to the ancestor of T. terrestris and in the western lowlands 
(Choco refugium) to the ancestor of T. bairdii, which later migrated 
to Central America (one unique lineage hypothesis). An alterna¬ 
tive hypothesis suggests that the ancestor of T. bairdii originated 
in a Central American refugium and then migrated towards the 
Pacific area of what would eventually become the western section 
of current Colombia (Haffer 1970). 

Ashley et al. (1996) used mitochondrial COII sequences to deter¬ 
mine that the ancestor of T. bairdii diverged from the ancestor of 
the two strictly South American tapir species around 19-20 MYA. 
The temporal split between the ancestors of T. terrestris and T. pin¬ 
chaque was estimated around 3 MYA. Norman and Ashley (2000) 
estimated new temporal splits. For the ancestors of T. bairdii and 
T. terrestris + T. pinchaque, they ranged from 18 to 20.4MYA (COII), 
or from 15 to 16.5 MYA (32S rRNA). These genes supported diver¬ 
gence times between the ancestors of T. terrestris and T. pinchaque 
of 2.5-2.7 MYA and 1.5-1.6 MYA, respectively. Ruiz-Garcia et al. 
(2012) used mitochondrial Cyt-b sequences and a Bayesian tree 
to determine that the ancestors of T. terrestris and T. pinchaque 
diverged around 3.8 MYA. Using a Median Joining Network, the 
most frequent T. terrestris haplotypes diverged from the main T. 
pinchaque's haplotype around 1.55 ± 0.32 MYA. 

Recently, Cozzuol et al. (2013) claimed the existence of a new 
tapir species. They compared Cyt-b gene sequences of four Amazon 
tapir individuals (two from Brazil and two from Colombia) with 
45 Cyt-b gene sequenced from Thoisy et al. (2010). Additionally, 
they included results for two mitochondrial genes (COI and COII) 
for six T. terrestris, one T. pinchaque and three individuals of T. kabo- 
mani. They concluded that the three mitochondrial genes showed 
the existence of T. kabomani as a full species. Voss et al. (2014) how¬ 
ever, questioned the validity of T. kabomani as a different species 
from T. terrestris. 

We sequenced the mitochondrial genomes of 45 T. pinchaque 
individuals sampled throughout six populations of Colombia and 
Ecuador. These two countries harbour more than 98% of the living 
mountain tapirs. 

The main aims of the present work are as follows: (1) To 
determine mitochondrial diversity levels in T. pinchaque; (2) To 
determine levels of genetic differentiation between T. pinchaque 


populations to provide new results testing for the possible exist¬ 
ence of subspecies within mountain tapirs (sensu Hershkovitz, 
1954); (3) To analyze possible historical demographic changes 
throughout the evolution of T. pinchaque ; (4) To determine spa¬ 
tial genetic structure in T. pinchaque; and (5) To contribute with 
new data to our understanding of speciation in the Tapirus genus. 

Material and methods 

Out of 80 T. pinchaque samples from Colombia and Ecuador we 
selected 45 (high-quality DNA) for mitochondrial DNA sequencing. 
Additionally, 17 samples of T. terrestris, 13 samples of T. bairdii and 
7 samples of T. indicus were also sequenced. The origins of sam¬ 
ples are in Table 1 and Fig. 1. T. pinchaque samples were obtained 
from hunted and captive animals as well as from tapirs captured 
for research (radio tracking). Some were from protected areas in 
Colombia and Ecuador (Table 2). 

DNA of high quality was extracted and isolated from blood and 
muscle samples using the QIAamp DNA Mini Kit (Qiagen, Inc.). 
Amplifications were carried out using the LongRange PCR Kit (Qia¬ 
gen, Inc.), with a final volume reaction of 25 pi. The reaction mix 
was composed of a 80-200 ng DNA template, 2 units of Long-Range 
PCR Enzyme, 3 pi of lOx LongRange PCR Buffer, 4 pi (15 pmol) of 
each primer, and 2 pi of 10 mM dNTPs. The cycling conditions were 
95 °C for 3 min, followed by 50 cycles denaturing at 95 °C for 20 s, 
primer annealing at 53-58 °C (with a decrease of 0.1 °C every cycle) 
for 30 s, and extension at 72 °C for 10 min. This was followed by 30 
cycles of denaturing at 95 °C for 20 s, annealing at 48-53 D C for 30 s, 
and extension at 72 °C for 5 min, with a final extension at 72 °C for 
10 min. All amplifications, including positive and negative controls, 
were checked onto 2% agarose gels under a Hoefer UV Transillumi¬ 
nator. Both mtDNA strands were sequenced directly using BigDye 
Terminator v3.1 (Applied Biosystems, Inc.). We used a 377A (AB1) 
automated DNA sequencer. 

We used four sets of primers (Table 3) to generate overlapping 
amplicons from 3,345 bp to 5,049 bp in length. These amplicons 
allowed us to carry out a quality test for genome circularity 
(Bensasson et al., 2001; Thalmann et al., 2004). Herein, we show the 
results of 15 mitochondrial genes (two rRNA and 13 protein codi¬ 
fying genes; Table 4). The sequences were concatenated by means 
of the SequenceMatrix v. 1.7.6 (Vaidya et al., 2011). They covered 
13,939 bp which represented about 83.1% of the total mitochondrial 
DNA length. Overlapping regions were examined for irregularities, 
such as frameshift mutations and premature stop codons. All the 
population genetics analyses carried out were using the sequences 
of these 15 mitochondrial genes. 

Additionally, we analyzed the phylogenetic relationships 
between T. pinchaque and T. kabomani using three concatenated 
mitochondrial genes (Cyt-b + COI + COII), the same genes used by 
Cozzuol et al. (2013). They published a tree including one sample 
of T. pinchaque, six of T. terrestris and three of T. kabomani. As T. kabo¬ 
mani was not included in the clade of T. terrestris + T. pinchaque, they 
concluded T. kabomani to be a new full species. To investigate this 
question, we sequenced 93 Neotropical tapir specimens (46 T. pin¬ 
chaque, 29 T. terrestris, 5 T. kabomani - including the two Brazilian 
specimens of Cozzuol et al., 2013, and 13 T. bairdii) for Cyt-b, COI 
and COII. The origins of the five T. kabomani are in Table 1. 

Population genetics analyses 

The Modeltest 3.7 (Posada and Crandall, 1998) and the Mega 5.1 
Softwares (Tamura et al., 2011) were applied to determine the best 
evolutionary nucleotide model for sequences of T. pinchaque, using 
the Akaike Information Criteria (A1C; Akaike, 1974) to determine 
the best evolutionary nucleotide substitution model. 
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Table 1 

Countries, number of individuals, exact geographic origins and the mitogenomes determined for the 45 Tapirus pinchaque, 17 Tapirus terrestris, eight Tapirus bairdii and seven 
Tapirus indicus individuals sequenced at 15 mitochondrial genes. 


Country 

Geographic origins 

Number of animals 

Mitogenomes found in each locality 

Tapirus pinchaque 

Colombia 

Los Nevados National Park in the Caldas Department 

3 

Ml, M2 


Los Nevados National Park in the Otun area in Risaralda Department 

3 

Ml, M2 


Los Nevados National Park in the El Bosque area in Risaralda Department 

2 

M2, M4 


Gaitania in Tolirna Department 

3 

Ml, M2.M15 


La Planada (La Bella and San Miguel), Marquetalia in Tolirna Department 

3 

M1.M19, M21 


Hereje River in Tolirna Department 

1 

M15 


Saldana River in Tolirna Department 

1 

M16 


La Azulena in Tolirna Department 

1 

Ml 


Penas Blancas in Tolirna Department 

1 

Ml 


Nevado del Huila National Park 

1 

M20 


Purace National Park in Huila Department 

1 

M3 

Ecuador 

Cuyuja, Quijos, Napo Province 

6 

M1.M6.M7, M12.M13, M14 


Oyacachi, Chaco, Napo Province 

3 

M1.M10, Mil 


Alto Papallacta River, Cayambe-Coca National Park, Napo Province 

9 

Ml, M5.M6, M7, M8, M9 


Cosangua, Napo Province 

2 

M13, M18 


La Bonita, Sucumbios Province 

1 

M3 


Sangay National Park (Las Culebrillas), Morona-Santiago Province 

3 

M12, M14.M17 


Podocarpus National Park, Loja Province 

1 

M13 

Tapirus terrestris 

Colombia 

Yali, Caqueta 

1 

M22 


Puerto Berrio, Antioquia 

1 

M23 


Remedios, Antioquia 

1 

M24 


Sabanacules, Antioquia 

1 

M25 

Ecuador 

Reserva Natural Cuyabeno 

1 

M26 

Peru 

Amaruyama, Madre de Dios River 

1 

M27 


Reserva Ecologica Tariyaca, Madre de Dios River 

1 

M28 


Parque Nacional Bahuaja-Sone 

1 

M29 


Tournavista, Pachitea River 

1 

M30 

Bolivia 

Ibare River 

1 

M31 


San Jose de Chiquitos 

1 

M32 


Trampa del Tigre 

1 

M33 

Argentina 

Chaco 

1 

M34 

Brazil 

Javari River 

1 

M35 

Suriname 

Magali 

3 

M36 

Tapirus bairdii 

Colombia 

Choco 

1 

M37 

Panama 

Darien 

8 

M37, M38, M39, M40, M41, M42 

Costa Rica 

Parque Nacional Braulio Carrillo 

3 

M43 

Mexico 

Yucatan 

1 

M44 

Tapirus indicus 

Malaysia 


3 

M45, M46 

Thailand 


3 

M46, M47 

Sumatra 


1 

M48 

Tapirus kabomani 

Brazil 

Southern Amazon Department 

1 

Cozzuol et al. (2013) 


Rondonia Department 

1 

Cozzuol et al. (2013) 

Colombia 

San Martin de Amacayacu, Amazon River 

1 

This work 

Peru 

Mazan River, affluent Napo River, Loreto, Northern Peruvian Amazon 

1 

This work 

Ecuador 

Tena, Upper Napo River, Ecuadorian Amazon 

1 

This work 


The following statistics to determine the genetic diversity for 
the overall sample of T. pinchaque and for the samples separately 
from Colombia and Ecuador were employed: number of haplotypes 
(H), haplotype diversity (H<j), nucleotide diversity (tt), average 
number of nucleotide differences (k) and 0 statistic by sequence. 
These genetic diversity statistics were calculated with the Programs 
DNAsp 5.1 (Librado and Rozas, 2009) and Arlequin 3.5.1.2 (Excoffier 
and Lischer, 2010). 

We used different tests to measure genetic heterogeneity and 
possible indirect gene flow estimates among the T. pinchaque popu¬ 
lations. Six geographical populations were considered, three in 
Colombia and three in Ecuador because there were six discrete 


geographical areas separated by considerable geographical dis¬ 
tances. We estimated the Gsx, 7st. N$x and Fsx statistics (Hudson 
et al., 1992). Significance was estimated with permutation tests 
using 10,000 replicates. 

We carried out an AMOVA analysis of the T. pinchaque sam¬ 
ple to determine the distribution of genetic diversity at different 
hierarchical geographical levels (Excoffier et al., 1992). Two groups 
were considered (Colombia and Ecuador) and within these groups, 
six populations (three in Colombia: Caldas-Risaralda, Tolirna and 
Huila; three in Ecuador: Napo, Sangay and Podocarpus). The fix¬ 
ation indices of Wright (1951) were estimated in the AMOVA 
analysis: <t> S c (variation of populations within the groups), <J> c t 
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Fig. 1. Map with the areas where 45 Tapirus pinchaque individuals were sampled in Colombia and Ecuador. 


(variation among groups) and <J> s t (variation among populations). 
Analyses were carried out using Arlequin 3.5.1.2 Software. 

A Bayesian skyline plot (BSP) was obtained by means of the 
BEAST v. 1.6.2 (Drummond and Rambaut, 2007) and Tracer vl.5 
Softwares to determine possible demographic changes across 
the natural history of T. pinchaque. The Coalescent-Bayesian 


skyline option in the tree priors was selected with three steps 
and a piecewise-constant skyline model with 60,000,000 gener¬ 
ations. Marginal densities of temporal splits were analyzed and 
the Bayesian Skyline reconstruction option was selected for the 
tree log files. A stepwise (constant) Bayesian skyline variant was 
selected with the maximum time as the upper 95% high posterior 


Table 2 

Protected areas (in km 2 ) in Colombia, Ecuador and Peru with confirmed presence of Tapirus pinchaque: estimations of low and high number of mountain tapirs in these areas, 
following Cavelier et al. (2010) and which of these protected areas were represented in our mitogenomic analysis. NP = National Park. 


Country 

Place 

Area (km 2 ) 

Lower estimation 

Upper estimation 

Analyzed in this study 

Colombia 

Los Nevados NP 

98,503 

17 

18 

Yes 


Las Hermosas NP 

301,045 

51 

55 

No, but some samples were 
from very nearby areas 


Nevado del Huila NP 

1,084,443 

185 

197 

Yes 


Purace NP 

618,197 

105 

112 

Yes 


Cueva de Los Guacharos NP 

39,816 

7 

7 

No 


Alto Fragua-Indi Wasi 

64,463 

11 

12 

No 


Sumapaz NP 

714,590 

122 

130 

No 


Cordillera de los Picachos NP 

264,495 

45 

48 

No 

Ecuador 

Cayambe-Coca NP 

1,408,606 

245 

261 

Yes 


Sumaco Napo Galeras NP 

525,865 

90 

95 

No 


Antisana NP 

617,042 

105 

112 

No 


Llanganatis NP 

824,778 

141 

150 

No 


Sangay NP 

1,294,180 

220 

235 

Yes 


Podocarpus NP 

1,069,749 

182 

194 

Yes 

Peru 

Tabaconas-Namballe 

238,907 

41 

43 

No 

Total 


9,194,679 

1,567 

1,669 
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Table 3 

Long range amplification primers employed in this work. 


Primer name 

Primer sequence (5' to 3') 

Amplicon sizes (bp) 

TPF1 

S'-CTAACTCCATAAAACA-S' 

3,768 

TPR1 

5'-CAAAATCCTCCGT-3' 


TPF2 

5'-AAACATAAGAAAT-3' 

4,990 

TPR2 

5'-ATCAGCCCTCCTTA-3' 


TPF3 

5'-TAGCCTTAATTCAA-3' 

5,049 

TPR3 

5'-TAAnTCCCCATTA-3' 


TPF4 

5'-TCAAACCAAAAAGG-3' 

3,345 

TPR4 

5'-CAGGTTTGATCCT-3' 



density (HPD). The time during which the demographic changes 
were recorded was chosen by the software with a mutation rate of 
2.5% per million of year (Ruiz-Garcia et al., 2015). 

We used IBD Software (version 1.2) (Bohonak, 2002) to deter¬ 
mine possible correlation between the genetic distances (Kimura 
2P distance; Kimura, 1980) and the geographical distances among 
these populations. A Mantel’s test (Mantel, 1967) was used to deter¬ 
mine the significance between the genetic distance and geographic 
distance matrices. The intercept and the slope of this relation¬ 
ship were calculated using Reduced Major Axis (RMA) regression 
(Sokal and Rohlf, 1981 ). Ten thousand randomizations were carried 
out to determine 99% confidence intervals. The calculations were 
completed with non-transformed and In transformed data (genetic 
distance & geographical distance). 

Phylogenetic analyses 

Two phylogenetic trees were obtained for the complete 15 
mitochondrial genes concatenated (13,939 bp). The first one 
was calculated with the maximum likelihood procedure (ML; 
Felsenstein, 1981) with the General Time Reversible model 
(GTR + G + I; G = gamma distributed rate variation among sites; 
I = constant proportion of invariable sites), being the best nucleotide 
substitution model in our data found by the Modeltest 3.7 Software 
(Ln = -25,474.80). This tree was obtained with the PAUP*4.0b8 Pro¬ 
gram (Swofford, 2002). The second one was a Bayesian tree (BI; 
Mau, 1996; Mau et al., 1999; Rannala and Yang, 1996). It was also 
based on the GTR + G + I model. This Bayesian analysis was carried 
out with the BEAST (vl.6.2) Software (Drummond and Rambaut, 
2007). 

Two runs of analyses were done, assuming a Yule speciation 
model and a relaxed molecular clock with an uncorrelated log¬ 
normal rate of distribution (Drummond et al., 2006). Results from 
the two independent runs (60,000,000 generations with the first 


Table 4 

Fifteen mitochondrial genes sequenced and sizes in base pairs (bp) for 45 Tapirus 
pinchaque, 17 Tapirus terrestris, eight Tapirus bairdii and seven Tapirus indicus. 


Mitochondrial genes sequenced 

Sizes (bp) 

I2S rRNA 

949 

I6S rRNA 

1,578 

NDI 

956 

ND2 

1,041 

ND3 

345 

ND4L 

296 

ND4 

1,377 

ND5 

1,820 

ND6 

526 

Cyt-b 

1,139 

COI 

1,544 

con 

683 

com 

783 

ATP8 

203 

ATP6 

680 


6,000,000 discarded as burn in and parameter values sampled 
every 10,000 generations) were combined with LogCombiner Soft¬ 
ware (vl.6.2) (Drummond and Rambaut, 2007). The final tree was 
estimated using TreeAnnotator Software (vl.6.2) (Drummond and 
Rambaut, 2007) and visualized in the FigTree vl.3.1 Program. Addi¬ 
tionally, the BEAST (vl.6.2) Software was run to estimate the time 
to the most recent common ancestor (TMRCA) for the different hap- 
lotype lineages found in different tapir species. As priors, we took 
a temporal split of 15 ± 2 MYA between the ancestor of T. indicus 
and those of the Latin American tapirs. We also used the temporal 
split of 7.5 ± 1 MYA between the ancestor of T. bairdii and those of 
the South American tapir species following Cozzuol et al. (2013). 
These estimates were recorded from the fossil record (Ferrero and 
Noriega 2007; Holanda et al., 2011; Holanda and Ferrero, 2013). 

Additionally, we created a maximum likelihood tree, for the 
93 tapirs sequenced at the concatenated Cyt-b + COI + COII genes to 
analyze the relationships among T. pinchaque, T. terrestris and T. 
kabomani. 

Results 

Genetic diversity 

The maximum likelihood estimate of Transition/Transversion 
was 2.37 (Ln = -26,308.12). For a GTR model, the maximum likeli¬ 
hood estimate of the gamma parameter for sites rates was 0.2037 
(+ G; five categories) (Ln =-25,481.66). 

The genetic diversity level for the global T. pinchaque sample 
analyzed for the 15 mitochondrial genes concatenated showed 21 
different haplotypes with H d =0.904± 0.003 (high genetic diver¬ 
sity) and tt = 0.0041 ±0.0003 (low-medium genetic diversity) (see 
Table 5). The genetic diversity levels for the Colombian and the 
Ecuadorian mountain tapir populations were similar, although the 
most unbiased gene diversity statistics were somewhat higher for 
the Colombian population. 

Genetic heterogeneity tests (Table 6A) revealed a low genetic 
differentiation level (G ST = 0.0231, 7 ST = 0.0534, N ST = 0.0558 and 
Fst = 0.056; P = 0.714-0.367, no significant), with the gene flow esti¬ 
mates showing that the populations were not disconnected and/or 
their separation was very recent (Nm = 1.58-9.15). F ST analysis 
(Table 6B) showed that only three population comparisons out of 
15 were significant. A gene flow comparison of population pairs 
only revealed 4 out of 15 cases with values lower than one. 

The AMOVA revealed that the percentage of variation between 
groups (Colombia and Ecuador) was very small (2.27%) and not 
significant (4>cr = 0.0227, P = 0.1026 ±0.0029). The percentage of 
variation among the populations within the groups was also small 
(6.76%) and not significant (4>sc = 0.0692, P = 0.1507 ±0.0038). The 
major percentage of genetic variance within populations was 
90.97%, but the difference between all the populations was not 
significant (4>st = 0.0903, P = 0.1381 ±0.0033). 

Historical demographic changes and spatial structure 

The BSP analyses (Fig. 2) supported a continuous population 
expansion for the total population throughout the last 100,000 
years, with the highest increase of females in the last 25,000 years. 
The Ecuadorian sample also yielded a continuous increase in female 
numbers within the last 25,000 years but this increase was more 
noteworthy in the last 10,000 years. The Colombian sample pre¬ 
sented a more constant population size in the last 100,000 years, 
with a strong decrease in the last 5,000 years. 

The four regression models did not show any significant rela¬ 
tionship between genetic and geographic distances. The Mantel test 
between the ln(genetic distance) and the ln(geographic distance), 
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Table 5 

Genetic diversity statistics for Tapirus pinchaque for the overall sample, for the Colombian sample and for the Ecuadorian sample at the 15 mitochondrial genes sequenced. 
Number of haplotypes (NH), haplotype diversity (H d ), nucleotide diversity (it), average number of nucleotide differences (I<) and 0 statistic (=2N e p,; N,. = effective female 
population size; p. = mutation rate per generation) by sequence. 



NH 

Hd 

TT 

I< 

0 per sequence 

Total sample studied of Tapirus pinchaque 

21 

0.904 ± 0.003 

0.0041 ± 0.0003 

4.457 ± 0.238 

15.094 ± 0.582 

Colombian sample studied of Tapirus pinchaque 

9 

0.816 ± 0.013 

0.0053 ± 0.0009 

5.726 ± 0.862 

13.248 ± 1.799 

Ecuadorian sample studied of Tapirus pinchaque 

14 

0.937 ± 0.013 

0.0029 ± 0.0007 

3.200 ± 0.710 

6.356 ± 1.366 


for instance, showed r =-0.127 (p< 0.691 from 10,000 randomiza¬ 
tions). 

Phylogenetics of T. pinchaque 

The ML tree (Fig. 3A) did not yield comprehensive clades reflect¬ 
ing geographical distribution, which agreed quite well with the 1BD 
analysis. There were only some small clades with animals of related 
geographical areas. These clades comprised five Colombian animals 
from Los Nevados National Park and Tolima (bootstrap 85%) and 
tree animals from the Ecuadorian Napo Province (49%). One individ¬ 
ual from La Planada (Tolima, Colombia) was the most differentiated 
of all the T. pinchaque studied. The B1 tree used T. indicus as outgroup 
and it showed the first split by the ancestor of T. bairdii (p = 1) and 
later the ancestors of T. terrestris and T. pinchaque (p = 1) (Fig. 3B). 
The two most divergent T. pinchaque were two individuals from La 
Planada-Tolima (Colombia) and Chaco-Oyacachi-Napo (Ecuador) 
(p = 0.9 and 1, respectively). The well-supported clades were the 
same as those revealed in the previous tree. Therefore, the phylo¬ 
genetic analyses revealed a low degree of geographic structure in 
the mountain tapir. 

Temporal splits within the genus Tapirus 

With the Bayesian procedure we determined a temporal sep¬ 
aration of the ancestor of T. bairdii relative to the ancestor of the 
South American tapirs ranging from 10.5 to 4.63 MYA (95% HPD; 
p, = 8.1 MYA). The temporal split between the ancestors of T. ter¬ 
restris and T. pinchaque oscillated from 7.42 to 3.27 MYA (95% HPD; 
p, = 3.7 MYA). In our BI tree, the mitochondrial haplotype diver¬ 
sification within T. terrestris began somewhere between 5.34 to 
2.36 MYA (95% HPD; p = 3.3 MYA). Within T. pinchaque, this process 
began around 6.09-2.68 MYA (95% HPD), with p = 2.9 MYA. 

Our ML tree with 93 specimens showed reciprocal monophyly 
between T. terrestris and T. pinchaque. (Fig. 4) Also, the mitochon¬ 
drial data support T. kabomani as a clade within T. terrestris. 


Discussion 

Genetic diversity, heterogeneity, demographic changes and 
absence of spatial structure in T. pinchaque 

Our nucleotide diversity values for T. pinchaque 
(tt = 0.004±0.0003) are considerably lower than estimates 
for T. terrestris reported by other studies. Ruiz-Garcia et al. (2015) 
estimated a value of tt = 0.0114 ±0.0003 for a wide geographical 
sample of T. terrestris. 

The overall sample of T. pinchaque has about three times lower 
genetic diversity than T. terrestris for mitochondrial genes. How¬ 
ever, T. pinchaque presented a higher nucleotide diversity than T. 
bairdii (tt = 0.0025 ±0.0005; Ruiz-Garcia et al., 2012), although the 
geographical distribution and population size of T. pinchaque are 
considerably smaller than in T. bairdii. Cavelier et al. (2010) esti¬ 
mated a maximum population of 5,000 to 5,700 mountain tapirs, 
based on existing suitable habitat, and estimated densities for this 
tapir species. However, the same authors considered that the pop¬ 
ulation censuses are really smaller than 5,000 animals, and that the 
overall mountain tapir population could be around 2,650 to 2,850 
individuals. 

Currently, seven (Lizcano et al., 2002) or eight (Cavelier et al., 

2010) protected areas contained mountain tapirs in Colombia 
and six protected areas contained mountain tapirs in Ecuador. 
Our molecular results, therefore, adequately represent the pop¬ 
ulation of T. pinchaque within the Colombian Central Cordillera. 
However, we could not sequence any individual from the East¬ 
ern Cordillera nor from the Western Cordillera (in the supposed 
case that this species lives in that area). It would be highly benefi¬ 
cial to analyze individuals of these Colombian cordilleras because 
their populations are probably extremely small and at an elevated 
risk of extinction. However, our results are representative of all 
the Ecuadorian T. pinchaque populations because we analyzed sam¬ 
ples from the three main areas where this species lives in Ecuador 
(Northern Eastern Cordillera: Cayambe-Coca National Park; 


Table 6 

Genetic heterogeneity and gene flow (Nm) statistics for different Tapirus pinchaque sets at 15 mitochondrial genes: (A) Six T. pinchaque populations (three in Colombian 
and three in Ecuador): (B) Population Fst (below) and Nm (above) pairs among the six T. pinchaque populations analyzed. *Significant Probability (P<0.05). inf = infinite. 
1 = Los Nevados National Park (Colombia); 2=Tolima Department (Colombia): 3 = Purace National Park (Colombia); 4 = Cayambe-Coca National Park and other areas in Napo 
Province (Ecuador); 5 = Sangay National Park (Ecuador); Podocarpus National Park (Ecuador). 

(A) Six T. pinchaque populations (three in Colombian and three in Ecuador) 

Statistics Gene flow 


Gst = 0.0518 

7st = 0.1612* 

Nst = 0.2389 ! 
Fst = 0.2401* 


Nm = 9.15 
Nm = 2.60 
Nm = 1.59 
Nm = 1.58 


(B) Population Fst (below) and Nm (above) pairs among the six T. pinchaque populations analyzed 



1 

2 

3 

4 

5 

6 

1 

0.000 

16.42 

0.34 

2.85 

0.79 

1.35 

2 

0.029 

0.000 

76.55 

8.80 

67.21 

inf 

3 

0.595* 

0.006 

0.000 

1.06 

0.64 

0.20 

4 

0.149* 

0.054 

0.319 

0.000 

4.75 

inf 

5 

0.387* 

0.007 

0.438 

0.095 

0.000 

inf 

6 

0.270 

0.000 

0.115 

0.000 

0.000 

0.000 
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Fig. 2. Bayesian skyline plot (BSP) analyses applied to the Tapiruspinchaque population studied at the mitochondrial DNA. (A) Overall Tapiruspinchaque sample; (B) Colombian 
mountain tapir sample; (C) Ecuadorian mountain tapir sample. On the x-axis, time in millions of years; on the y-axis, size of female effective population. Dashed vertical 
lines are initial points where demographic changes were detected. 
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Fig. 3. (A) Maximum likelihood tree with 45 Tapirus pinchaque, plus 17 7. terrestris by employing sequences of 15 concatenated mitochondrial genes. (B) The same with a 
Bayesian tree. The number in the nodes are bootstrap percentages higher than 50% for the first tree; the number in the nodes are posterior probabilities higher than 0.5 for 
the second tree. 
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Fig. 3. (Continued) 


Central Eastern Cordillera: Sangay National Park; and Southern 
Eastern Cordillera: Podocarpus National Park). Thus, our genetic 
diversity estimates for the mountain tapir population of Ecuador 
should be more accurate and realistic than what we estimated for 
the Colombian population. 


Our BSP analysis indicated a strong population decrease in the 
last 5,000 years for the Colombian population, which coincided 
with one of the dry periods in the Holocene after the Optimum Cli- 
maticum. Following Rothlisberger (1987), around 6.300YA, there 
was a significant increase in temperature especially in Southern 
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T. pinchaque 


T. bairdii 


Fig. 4. Maximum likelihood tree for 93 Neotropical tapir specimens (including T. pinchaque, T. terrestris, T. bairdii and T. kabomani ) for three concatenated mitochondrial 
genes (Cyt-b + COI+COII). 


South America as well as in the Central Andes. This dry period, 
around 5,000-6,000 YA, was also detected in the Amazon, Caqueta 
and lower Magdalena River basins as well as in the Andean lagoons 
of Colombia and Peru (Thompson et al„ 1995; Van der Hammen, 


2001). These climatic changes together with hunting by humans 
(Frisson, 1998) might explain the population decline for T. pin¬ 
chaque in the area of Colombia. This fact also coincides with that 
found by Ruiz-Garcia et al. (2015) for the T. terrestris population of 
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Northern Colombia and Northwestern Venezuela, which showed 
an extreme population decrease in the last thousand years. In 
contrast, the Ecuadorian population showed a relatively recent 
population expansion (around 10,000 YA). The most intense cold 
and dry period of the last Pleistocene glaciation occurred 26,000 
to 14,000 YA (Van der Hammen and Cleff, 1992). This period was 
related to the maximum extension of the glaciers in the north¬ 
ern hemisphere and the complete disappearance of some lagoons 
near Bogota (Van Geel and van der Hammen, 1973). Later, from 
14,000 YA to 10,000 YA (Tardiglacial), the climate became hotter 
and wetter. After these dry and cold periods the mountain tapirs 
in Ecuador probably experienced a population increase due to bet¬ 
ter climatic conditions. It seems that the climatic Holocene event 
which affected the mountain tapirs in Colombia did not affect the 
tapirs in Ecuador. An unequal impact of hunting in these areas could 
also explain these differences. 

Non-significant or very limited genetic heterogeneity among 
mountain tapir populations means that historical gene flow has 
been sufficient enough to avoid genetic differentiation, although 
the Andean habitat is currently extremely fragmented. Fragmenta¬ 
tion might be very recent, thus preventing genetic differentiation 
of the mountain tapir populations. Another possibility might be 
high gene flow. This has been demonstrated for T. terrestris (Pinho 
et al., 2013 ) but not yet for T. pinchaque. The non-existence of spatial 
genetic structure is an interesting fact from a conservation perspec¬ 
tive because there should not be genetic consequences if animals 
are translocated. 

Phylogenetics of the Neotropical tapirs and evolutionary patterns 
in T. pinchaque 

Hershkovitz (1954), Haffer (1970) and Holanda and Ferrero 
(2013) used morphological characters and determined that the lin¬ 
eage which gave rise to T. pinchaque was the oldest of the current 
Neotropical tapir species. The conclusion by Hershkovitz (1954) 
was based on the close similarity of cranial contours to some fossil 
tapir skull which has a less specialized proboscis than current tapir 
species. 

The adaptation of T. pinchaque to the temperate and cold cli¬ 
mates motivated Hershkovitz (1954) to affirm that this species 
appeared in a temperate climate before the completion of the Pana¬ 
manian land bridge with the minimum age of this species being 
around 3-4 MYA. Lack of fossil records of T. pinchaque makes it dif¬ 
ficult to support this hypothesis. Additionally, no molecular results 
support the fact that the ancestor of T. pinchaque is the most ancient 
of the current Neotropical tapirs. All the molecular studies have 
shown that T. pinchaque and T. terrestris are sister groups (Ashley 
et al„ 1996; Norman and Ashley, 2000; Ruiz-Garcfa et al., 2012, 
2015). 

The problem of using craniometrical and teeth data when deal¬ 
ing with recently differentiated organisms is that nobody knows 
the degree of phenotypic plasticity of these characters and their 
genetic basis. For instance, Hershkovitz (1954) found a variable 
trait in the first upper premolar of T. pinchaque. Ecuadorian speci¬ 
mens show the simple condition, with the cinguloid shelf absent, 
whereas the first Colombian skull analyzed of T. pinchaque had a 
premolar like another Tapirus species. For this reason, two sub¬ 
species (T. p. pinchaque and T. p. leucogenys ) were identified within 
the mountain tapir. In contrast, our mitogenomic analysis shows 
that haplotypes from animals from both distributions are inter¬ 
mixed and it could represent a unique taxon and/or they diverged 
very recently. Molecular data do not support the paleontological 
hypothesis of Holanda and Ferrero (2013) nor the hypothesis of 
Hershkovitz (1954) that T. pinchaque and T. terrestris correspond to 
different lineages of tapirs. Our mitogenomic results are more in 
line with the paleontological view of Hulbert and Wallace (2005) 


and Hulbert (2010) because they consider that T. pinchaque and T. 
terrestris form a monophyletic group. 

Another question is that the genetic distances between T. 
pinchaque and T. terrestris are very small for traditional fully differ¬ 
entiated species. Kartavtsev (2011) analyzed sequences of COI from 
20,731 vertebrate and invertebrate animal species and estimated 
the average distance data for five different groups. He calculated 
0.89%± 0.16% for populations within species and 3.78%± 1.18% for 
subspecies or semispecies. For species within a genus he calculated 
11.06% ± 0.53%. For this gene, the genetic differentiation between 
T. terrestris and T. pinchaque was only 1%, a value for populations 
or subspecies within species. However, we consider T. pinchaque a 
full species by different criteria. Using the PSC, it formed a recipro¬ 
cally monophyletic clade with T. terrestris (including T. kabomani ) 
and using the BSC, it is a reproductive unit because no natural 
nor captivity hybridization with T. terrestris has been reported due 
probably to important differences in their karyotypes, although its 
molecular differentiation is minor. Barongi (1993) and Houck et al. 
(2000) showed 2N = 80 and FN = 80-92 for T. terrestris, but 2N = 76 
and FN = 80-84 for T. pinchaque. The correlation of important chro¬ 
mosome differences but with very low mitochondrial divergence is 
typical of peripatric or parapatric chromosomal speciation (Lewis, 
1966; White, 1968, 1978; King, 1993). This has been reported in 
other Peryssodactyla, such as in the different zebra species (Bush, 
1981). Additionally, ecological isolation could be another cause of 
speciation for T. pinchaque (Ecological Species Concept). 

But when did the last common ancestor of T. terrestris and 
T. pinchaque live? Although, Cozzuol et al. (2013) (see Fig. 
7 of these authors) claimed that this could have happened 
around 0.1-0.3 MYA, we believe it occurred considerably before 
(1.5-7MYA) for two reasons: (1) Ashley et al. (1996) and Norman 
and Ashley (2000) estimated a split around 1.5-3 MYA with two 
mitochondrial genes. Ruiz-Garcla et al. (2012) showed a temporal 
split between the ancestors of both tapir species of around 3.8 MYA 
for a Bayesian tree with some temporal priors. With other proce¬ 
dures, without priors, these authors estimated a split of around 
1.6 MYA. Our present results with mitogenomics showed a range 
of 7.4-3.3MYA with a mean around 3.7 MYA. (2) Downer (2001) 
showed that a significant portion of the original high Andean flora 
from Northern Peru and further north could have depended upon 
the mountain tapir as a seed disperser. This Flora appeared during 
the geologically recent Andean rise that began after the Quaternary 
Period (late Cenozoic) around 2.5 million MYA (Simpson, 1979; 
Benton, 1991). This author showed indirect ecological evidence in 
favor of the presence of the mountain tapir in the last 2.5 MYA 
within the Northern Andean Cordilleras. 

Molecular evidence against T. kabomani as a new species 

Our results showed reciprocal monophyly between T. pinchaque 
and T. terrestris and that T. kabomani was a clade nested within 
T. terrestris. Thus, our mitochondrial analyses did not validate T. 
kabomani as a new species as suggested by Cozzuol et al. (2013). 
Our results agree with those of Voss et al. (2014), who claimed that 
T. kabomani is not a different species from T. terrestris. Cozzuol et al. 
(2013) analyzed a very small sample of T. pinchaque that probably 
did not represent the total mitochondrial genetic diversity of this 
species. By chance, it was nested within the genetic diversity of T. 
terrestris because of the small genetic distances between both taxa. 

Additionally, we did not detect a correlation between the mor- 
photype proposed by Cozzuol et al. (2013) for T. kabomani and 
the T. kabomani mitochondrial haplotypes. Cozzuol et al. (2013) 
described T. kabomani with a lower sagittal crest, a smaller size and 
a darker color than the typical morphotype of T. terrestris. How¬ 
ever, the Colombian (Amacayacu River) and the Ecuadorian (upper 
Napo River) animals which showed haplotypes of the T. kabomani 
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presented typical morphotypes of T. terrestris (we couldn’t deter¬ 
mine the morphotype of the Peruvian specimen because the sample 
came from a piece of skin of a hunted animal). Contrarily, we have 
found and sampled small-sized and dark lowland tapirs similar to 
the morphotype described by Cozzuol et al. (2013) for T. kabomani 
but they presented haplotypes of other T. terrestris haplogroups. 

We agree with Zachos (2015) who stated that species are 
such fundamental units that they should not be introduced care¬ 
lessly and that species description and splitting based on simple 
morphometric differences (even significant ones) or phylogenetic 
relationships derived from limited molecular datasets (for instance, 
only mtDNA) should be strongly discouraged. They may serve 
to support conclusions derived from larger and more complete 
datasets, but are not enough on their own. Therefore, additional 
data are needed to test if T. kabomani is a new species. Nuclear DNA 
should be analyzed to confirm the findings revealed by mitoge- 
nomics for T. pinchaque. 
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